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ABSTRACT 

Among efforts to detect gravitational radiation, pulsar timing arrays are uniquely poised to de¬ 
tect “memory” signatures, permanent perturbations in spacetime from highly energetic astrophysical 
events such as mergers of supermassive black hole binaries. The North American Nanohertz Obser¬ 
vatory for Gravitational Waves (NANOGrav) observes dozens of the most stable millisecond pulsars 
using the Arecibo and Green Bank radio telescopes in an effort to study, among other things, grav¬ 
itational wave memory. We herein present the results of a search for gravitational wave bursts with 
memory (BWMs) using the first five years of NANOGrav observations. We develop original methods 
for dramatically speeding up searches for BWM signals. In the directions of the sky where our sen¬ 
sitivity to BWMs is best, we would detect mergers of binaries with reduced masses of 10 9 M 0 out 
to distances of 30 Mpc; such massive mergers in the Virgo cluster would be marginally detectable. 

We find no evidence for BWMs. However, with our non-detection, we set upper limits on the rate at 
which BWMs of various amplitudes could have occurred during the time spanned by our data-e.g., 

BWMs with amplitudes greater than 10 -13 must occur at a rate less than 1.5 yr _1 . 

Subject headings: pulsars, gravitational waves 
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1. INTRODUCTION 

Due to the intrinsically nonlinear nature of Ein¬ 
stein’s equations, all systems that radiate gravitational 
waves (GWs) are anticipated to produce “memory,” non- 
oscillatory components of the gravitational waveforms 
( Sma rrl 19771 i Bontz fc Prlcep 979; Br aginskii fc Thorn e 
/987| |Christodoulou||l99l[ |Blancnet fc Damour||1992p. 
Supermassive black hole binaries (BlVlBHBs), during the 
final few orbits preceding their merger, are expected to 
generate GW bursts with memory (BWMs) with suffi¬ 
ciently large amplitudes to make them potentially de¬ 
tectable with pulsar timing arrays (PTAs; |Fav ata 2009; 
Seto||2009l IPshirkov et al.J 2010Uvan Haasteren fc'Levm 


2010[|Cordes fc Jenet|2012||Madison et al.|2014|) . |Cutler 

ently singled out memory as a key detec- 


et al 


tion target for PTAs as a probe of exotic and unexpected 
GW sources like phase transitions in the early Universe. 

To facilitate GW detection, several international con¬ 
sortia are currently using sensitive radio telescopes paired 
with pulsar-optimized hardware to realize the PTA con - 
cept (iHellings & Downs||1983| |Foster fc Ba cker|| 19901). 
The E nropean~Pnlsar Timing Array (EPTA; Kramer fc 
Champion|2013 ), the North American Nano h ertz ObseU - 
vat ory fo r Gravitational Waves (NANOGrav; McLaugh¬ 


lin 


2013), and the Parkes Pulsar Timing Array (PPTA; 

Mlt2( 


Hobbs 2013) are pushing precision pulsar timing to its 
limits and developing new data analysis techniques to 
usher in the era of PTA GW astronomy. By pooling their 
data and expertise, these consortia have formed_the In¬ 
ternational Pulsar Timing Array (IPTA; [Manchester & 
IPTA|2013|) which is poised to become the most sensitive 
of all PTAs. 
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In recent years, the various PTAs have begun to place 
astrophysically meaningful upper limits on conti nuous 
GWs from individually resolvab le SMBHBs (Arzouma- 
|nian et ah]|2014| Zhu et al. 20 141) and a stochastic back¬ 
ground of GWs (van Haaste ren et al. 2011; Dem orestl 
et al. 2013 Shannon et al. '2013k. Wang et al. (|2015|) 
have searched for BWIVIs in the first approximately six 
years of PPTA data; they detected nothing, but deter¬ 
mined with 95% confidence that BWMs with amplitudes 
greater than 10 -13 occur at a rate less than 0.8 yr -1 . 

In this paper, we search for BWMs in the first approx¬ 
imately five years of NANOGrav data using new tech¬ 
niques that lead to considerable computati onal speedups 
over t he search methods described by Madison et al. 
(|2014h and those used by [Wang et al.|| (20151). In Sec¬ 
tion 2, we describe the data used in our~BWM search. 
In Section 3, we describe the BWM signal model we use 
for our analysis. In Section 4, we discuss models of noise 
and how our sensitivity to BWMs is influenced by them. 
In Sections 5 and 6, we describe our search techniques, 
differentiating between searches for so-called pulsar-term 
and Earth-term events. In Sections 7 and 8, we present 
the results of our pulsar-term and Earth-term searches, 
respectively. In Section 9, we place upper bounds on 
BWM rates and amplitudes. In Section 10 we summa¬ 
rize our key results and offer concluding remarks. We 
provide a summary of our important notational elements 
in a table of key symbols in an appendix. 

2. PULSAR TIMING DATA SET 

In this section, we discuss several aspects of the data 
that are relevant to our a nalysis. Fo r a more t horough 
description of th e data, see Demorest (2007) and Demor- 
|est et al.| (2013). 

The five-year data set consists of pulse times-of-arrival 
(TOAs) collected at approximately monthly intervals for 
each of 17 millisecond pulsars (MSPs). All observations 
were done with the Arecibo radio telescope and the Green 
Bank Telescope (GBT). Of these 17 MSPs, those visi¬ 
ble to Arecibo were observed with Arecibo; all others 
were observed with the GBT. One pulsar, J1713+0747, 
was observed with both telescopes. All observations were 
done using one of two identical backend systems: the As¬ 
tronomical Signal Processor (ASP) and the Green Bank 
Astronomical Signal Processor (GASP). These backends 
performed real-time coherent dedispersion over bands up 
to 64 MHz wide and recorded the results averaged over 
channels of width 4 MHz each. 

For each observing epoch, several TOAs are reported 
from various frequency channels of the 64 MHz band. 
At Arecibo, observations were typically conducted at 
two widely-separated frequencies (usually 430 MHz and 
1400 MHz) within one day; at the GBT, observations at 
two frequencies (usually 820 MHz and 1400 MHz) were 
conducted within several days of each other. Approx¬ 
imately contemporaneous observations at multiple fre¬ 
quencies allow epoch-to-epoch timing fluctuations caused 
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by changes in dispersion m easu re (DM) to be accounted 


,_ o 

for in the timing model fit (Lam et al. 2014). 

Three of the pulsars comprising the five-year 
NANOGrav data set, J1853+1308, J1910+1256, and 
B1953+29, do not have sufficient dual-frequency cover¬ 
age to correct for timing errors from variations in DM 
over time; we exclude these pulsars from our analysis for 
this reason. We also exclude the data set for J1600—3053 
from our analysis because it is comparatively very short, 
spanning just two years. Searching for BWMs in such 
short data sets is feasible in principle, but we avoid it for 
two reasons. First, the minimum detectable BWM am¬ 
plitude in a particular data set scales app roximately as 
T~ 3 / 2 where T is the span of the data set (van Haasteren 
& Levin 2010), so we do not anticipate that this data set 
will greatly improve our sensitivity to BWMs. Second, 
in such short data sets, many timing model parameters 
are highly covariant with each o ther (e.g. , spin and as¬ 
trometric parameters; |Madison et al. 2013) and with any 
BWM signal present. 

Finally we exclude the data for J1643—1224 from our 
analysis because we believe it contains chromatic tim¬ 
ing biases from unaccounted-for phenomenology in the 
interstellar medium (ISM; see Figure 2). The DM of 
J1643—1224 is approximately 62 pc cm -3 , nearly a factor 
of two greater than any other pulsar in our sample. With 
high DM pulsars, chromatic timing errors that deviate 
from the u~ 2 scaling expected from cold p lasma dis per 


sion alone become more significant (Cordes & Shannon 
2010). Furthermore, this pulsar is directly behind a com¬ 
plex region of H-II associated with £-Ophiuchi, a mas¬ 
sive^ runaway O- type star spinning very near b reakup 
(Gaustad et al.|p 

H-Il 


2001 Villamariz & Herrero 


This 


intervening H-Il region could conceivably contribute to 
non-trivial and currently unaccounted-for chromatic ef¬ 
fects on the timing behavior of J1643—1224. 

3. SIGNAL MODEL 

For a given pulsar, the timing perturbation from a 
BWM of amplitude he is well-modeled as 

At(t) = h B B(6 ,0) [(t - t o )0(t - t 0 )- 

(t - h) 0(t - +)]. (1) 

The function B(0,(f)) = (1/2)(1 — cos 6) cos (20) ranges 
between —1 and 1 and is common to all pu lsar timin g 
efforts to detect point-like sources of GWs ( Estabrook 
fc Wahlquist||197'5l |Hellings fc Downs 1983|]Eee et al. 
2011). The angle between the direction the burst prop¬ 
agates and the line of sight from Earth to the pulsar is 
0; 4> is the angle between the principal polarization vec¬ 
tor of the wave and the projection of the line of sight 
from the Earth to the pulsar onto the plane normal to 
the wave propagation direction. The BWM encounters 
the Earth at a time to and is observed from Earth to 
encounter the pulsar at a time t\ = to + (Z/c) (1 + cos 0) 
where l is the distance from the Earth to the pulsar |van 


Haastere n fc Levin||2010| |Cordes fc Jenet||2012| [Madison 

et al.||2014p. The function 0 is the Heaviside step func- 

tion. The amplitude of a BWM coming from a SMBHB 
merger of reduced mass /i = M\M 2 /(Mi + M 2 ) ( M\ and 
M 2 are the masses of the black holes in the binary) with 
a typical inclination angle of X = n/3 at a luminosity 
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distance D B from Earth is 


( Madison et al.|| 2014) 


hs 


1.5 x 1(T 12 


io 9 M e J 


/1 Mpc\ 


( 2 ) 


The distances to the pulsars in the NANOGrav array 
are on the order of kiloparsecs. Our timing baseline is 
approximately five years. Unless 0 differs from 7r by less 
than ^ 3 degrees for a particular pulsar, we expect that 
to and t\ will not both fall within our observing window. 
Since we consider only 12 pulsars in our analysis, less 
than 1% of the sky is within 3 degrees of one of our pul¬ 
sars. Assuming BWMs occur isotropically, there is a less 
than 1% chance of both to and t\ occurring within our 
five-year observing span if a BWM occurs at all. And 
if a BWM occurred with a small enough angular sepa¬ 
ration from a pulsar in our array that we could see the 
timing perturbation both turn on at to and turn off at ti, 
it would only be observed to turn off in that one pulsar. 
So, in each of our pulsar timing data sets, we need only 
to look for evidence of timing perturbations of the form 


At(t) = h p (t - t B )0(t - t B ), (3) 

where h p = ±h B B(6, 0), what we call the projected burst 
amplitude, and t B is either to or ti, what we call the burst 
epoch. Bursts arriving at an individual pulsar only influ¬ 
ence the timing behavior of that pulsar; we refer to these 
as pulsar-term bursts (see Section 5). Bursts arriving 
at the Earth will simultaneously begin to influence the 
timing behavior of all pulsars in our sample; we refer to 
these as Earth-term bursts (see Section 6). 


4. NOISE MODEL 

Consider TOAs measured using pulse profiles obtained 
by synchronously averaging Af pulses for each of several 
radio-frequency channels. The pulse profiles are func¬ 
tions of pulse phase, ip, channel frequency, v , and ob¬ 
serving epoch, r. A TO A from a particular frequency 
channel and observing epoch can be written as 

tu,r = too,r + + tc u ,. r + 6 s/n ut + £j V jT + cdiss^j 

( 4 ) 


that the radiometer noise adds to a fixed pulse shape 
under the assumptions of matched filtering against a very 
high S/N template pulse profile. 

The second random contribution to t VjT , ej UT , is from 
phase jitter in single pulses. Jitter is highly correlated 
between frequency channels, but is known to decorre¬ 
late over widely-separated frequencies; the probability 
distribution function of the phase of single-pulse cen¬ 
troids evolves with fr equencies similar to pulse profiles 
( Shannon et al.||2014| ). Jitter noise can be modeled as 

( e Jv,T e J„r, T r ) = « Oj{N) 5 tt /, ( 6 ) 


where aj <x M 1 / 2 . The quantity pj t v , is approximately 
unity unless v and z/ are very widely separated. 

Usually, aj < (j g /n • Jitter noise begins to dominate 
radiometer noise only when the S/N of single pulses be¬ 
gins to exceed unity. For our data, collected with ASP 
and GASP, jitter noise should be significantly subdomi¬ 
nant to radiometer noise for all pulsars. Both radiometer 
and jitter noise are significantly larger than contributions 
from cdiss^ t , a random error associated with diffrac- 
tive interstellar sc intillation (DISS; |Cordes||1990| [Cordes 
fc Sha nnon 2010|). In a d etailed study of J1713+0747, 
Shannon & Cordes (2012) found that for single pulses, 


a j ~ 27±1 /is, while timing errors from DISS were only a 
few n anoseco nds. A recent 24-hour study of J1713+0747 
dDolch et al. |2014|) confirm ed the jitter measurement of 
Shannon & Cordes (2012). Errors from DISS can be 
substantially more important in the timing of high DM 
pulsars like B1937+21. We expect DISS to be a sub¬ 
dominant component of our timing error budget for the 
MSPs we include in our sample, so we ignore it in our 
analysis. In summary, the noise covariance matrix for 
a pulsar in our array, with these anticipated sources of 
noise, can be approximated as 


— ~ Jtt' ^S/N (AO T CFj(AT) 


( 7 ) 


where e v , T is the sum of all noise influencing the TOA 
from observing epoch r and frequency channel v. 


where t oo ?T is the TOA at infinite frequency, £dm is the 
dispersive delay from propagation through ionized inter¬ 
stellar plasma, tc is an additional, non-dispersive chro¬ 
matic perturbation, such as intrinsic profile evolution 
with frequency, pulse broadening from mu ltipat h prop 


agati on, and interstellar refraction (Cordes & Shannon 


2010). If unaccounted for, £dm and tc can lead to sys¬ 


tematic errors in TOA estimates. With multi-frequency 
observations at each observing epoch, effects from £dm 
are mitigated in the NANOGrav data set. Potential ef¬ 
fects from tc are combatted by fitting for constant inter¬ 
channel offsets unique to each pulsar. Unlike £dm and 
tc , the V terms in Equation 4 are random errors from a 
variety of noise sources. 

The term £s/n ut is uncorrelated between frequency 
channels and between observing epochs, i.e., 


(^/N^/S/N^y) — °S/n(A0 (5) 

where os/n oc ( S/N ) -1 oc A f -1 ' 2 . The quantity <ts/n 
is the TOA uncertainty from radiometer noise assuming 


4.1. Empirical Noise Models 


In a previous a nalysis of the first five years of 
NANOGrav data, Arzoumanian et al. (2014]) searched 
for continuous GWs using a different noise model: 


r' = 


Sw'(Q 2 + E 2 *£ /n (AT)) + J 1 


Cred (^5 T) • 


( 8 ) 


The terms Q and E are commonly referred to as EQUAD 
(a source of Gaussian white noise with time units added 
in quadrature to radiometer noise) and EFAC (a dimen¬ 
sionless constant multiplier on the anticipated amount 
of radiometer noise), respectively. The term J mimics 
the correlations induced between TOAs from different 
frequency channels of the same observing epoch by jit¬ 
ter, but qualitatively differs from jitter in that its value 
does not change to reflect the number of pulses, Af, aver¬ 
aged together to yield a TOA. Most observations in the 
five-year NANOGrav data set we considered had similar 
integration times, so Af does not fluctuate substantially 
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between epochs and J is thus very much like a jitter- 
induced correlation. Different values of Q, E , and J 
were used for each widely separated observing frequency 

of eac h pulsar. _ 

The |Arzoumanian et al. (2014) noise model also in¬ 
cluded red noise with a power-law power spectrum 
P(f) = A[f/(1 yr -1 )] 7 . This component of the noise 
model is not to be thought of as part of the TOA er¬ 
ror budget; TO As measured with extreme accuracy may 
still differ from the expectations of a timing model be¬ 
cause of spin noise intrinsic to a pulsar, which is well 
model ed as a st oc hastic pro cess with a p ower-law spec¬ 
trum dShan non fe Cordes 2010). In the Arzoumanian 
et al.| (|2014| analysis, under the assumption that no GW 
signal was present in the data set, the noise model pa¬ 
rameters 3 = [A, 7, E, Q, J] were determined by finding 
the maximum of the likelihood function 

= e*p[-l(R-Mip)rc-(R-Mip)] 

y (27r) JVTOA det C 

(9) 

where R are the timing residuals from an initial timing 
model, 5p are deviations of the timing model parameters 
from the initial timing model, M is the timing model 
design matrix, and TVtoa is the total number of TOAs 

in the data set. _ 

We applied t he noise model assessment done by!Arzou¬ 


manian et al. (2014|) to two simulations of the five-year 
NANOGrav data set for PSR J1909—3744, one that con¬ 
tained only simulated radiometer noise and one that con¬ 
tained radiometer noise and a bright BWM occurring at 
the midpoint of the data set. The resulting best-fit noise 
models were identical except for in the red noise parame¬ 
ters A and 7. The BWM signal is heavily covariant with 
the red noise component of the noise model, resulting in 
a large value for A and a poorly constrained value of 7. 
An indepe ndent analysis o f the fi ve-year NANOGrav 


data set by Perrodin et al. (2013) found that except 
for J1643—1224 and J1910+1256, pulsars we already ex¬ 
cluded from our analysis, all pulsars were consistent with 
white noise alone. Because the demonstrated covariance 


cate detection of a BWM ( 

Cordes & Jenet 

2012 

) and 

because of the lack of strong 

2; evidence tor rec 

1 noise, we 

have opted to consider the 
fixed noise models without t 

Arzoumanian et 

al. 

(2014) 

he red noise component. 


4.2. Comparing Noise Models 

In Figure 1, we show how the uncertainty on the pro¬ 
jected BWM amplitude, 07^, varies in the NANOGrav 
five-year data set for PSR J1713+0747 over several trial 
burst epochs under three different noise models. The 
quantity 07^ is a direct measure of how sensitive a par¬ 
ticular data set is to BWMs and it is dire ctly influ enced 
by the noise m odel; it is discussed bv, e.g., van Haasteren 
& Levin| (|2010|) and|Madison et al. (2014). 

The bottom curve of Figure 1 is based on a noise model 
that only includes radiometer noise. The middle curve 
uses the noise model described by Eq uati on 7 w ith the 
jitter measurement from Shannon & Cordes (2012). The 
top c urve i s based on the noise model used by |Arzo uma¬ 
nian et al.| (|2014|) and described in Equation 8 (without 
tfieTed noise component). The discrepancy between the 



Fig. 1.— Uncertainty on the projected amplitude of a BWM 
at various trial burst epochs with the NANOGrav five-year data 
set for PSR J1713+0747 using three different noise models. The 
lowest curve assumes that only radio meter n oise is p r esent in the 
data and is identical to the results of [Madison et al. | (2014). The 
middle curve assumes that only radiometer noise and jitter noise 
are present in the data (as in Equ ation 7); the scale of th e jitter 
contribution is consistent with |Shannon & Cordes (2012). The 
most conservati ve curve is b ased on the hxed noise model used 
in | Arzoumanian et al. | (|2014|) (as in Equation 8) sans a red noise 
component. 


physically motivated noise model of Equation 7 and the 
empirical noise model of Equation 8 is not currently well 
understood. Detailed studies like Dolch et al. (2014) are 
being carried out to better understand the noise budget 
of pulsar timing experiments and any systematic effects 
that influence the timing procedure, but bridging the gap 
between the top and middle curves of Figure 1 is an o n- 
going area of research. The Arzoumanian et al. (2014) 
noise model is the most conservative of the three we con¬ 
sider so we use it for the remainder of our analysis. 


4.3. Epoch Averaging 

In Figure 2, we depict the epoch-averaged timing resid¬ 
uals from the 12 pulsars in our sample (and the excluded 
pulsar J1643—1224). Residuals are obtained individually 
for multiple frequency channels for each integration last¬ 
ing 15 to 45 minutes. The residuals from that integration 
are combined to form a single epoch-averaged residual. 
The noise model is central to this procedure. We define 
an operator, 


A = (U T C _ 1 U ) _1 U T C _1 


( 10 ) 

which maps the raw residuals, R, to epoch-averaged 
residuals, R e = AR. The matrix U is the “exploder” 
matrix discussed in Arzoumanian et al. (2014) that maps 
epochs to the full set oFTOAs. The uncertainties on the 
epoch-averaged residuals are the square roots of the di¬ 
agonal entries of the matrix 


c E = (RfiR^) = (lEC- 1 !!) 1 . 


( 11 ) 


5. SEARCHES FOR BWMS IN THE PULSAR TERM 

Information in the pulsar terms of the 12 pulsars in our 
sample comes from causally distinct regions of spacetime 
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Fig. 2.— Epoch-averaged timing residuals for the 12 pulsars used 
in our analysis (and the excluded pulsar J1643—1224). The various 
colors indicate different observing frequencies: 327 MHz is green, 
430 MHz is blue, 820 MHz is red, 1400 MHz is black, and 2300 
MHz is cyan. We show the residuals for J1643—1224 to illustrate 
the apparent chromatic problems with the data: the residuals from 
820 MHz and 1400 MHz are consistently anti-correlated and are 
often extreme outliers. 


(|Cordes & Jenet|2012|. A BWM is 12 times more likely 
to encounter a single pulsar in our array than it is to 
encounter the Earth. However, it is not true that BWMs 
encountering individual pulsars in the array are 12 times 
more likely to be confidently detected. Different pulsars 
are timed with varying degrees of precision and with dif¬ 
fering observing cadences making them unequally sensi¬ 
tive in searches for BWMs. A particular BWM may be 
polarized in such a way or from such a part of the sky 
that the projection factor, £>(#,</>), makes its influence 
on the timing behavior of a particular pulsar vanishingly 
small. Furthermore, non-BWM pheno mena such as in¬ 
trinsic pulsar spin noise ( Shannon fe C ordes 2010) and 
microglitches (Cognard & Backer 2004| can be confused 
as BWMs. Without the signal appearing concurrently in 
multiple pulsars as in an Earth-term BWM, it is diffi¬ 


cult to rule out pulsar-specific phenomena. Nonetheless, 
we carry out a search for BWMs in each of our individ¬ 
ual pulsars as non-detections in many pulsars allow us 
to place constraints on otherwise inaccessible regions of 
BW M parameter space. 


Madison et al. (2014) describe techniques for search¬ 
ing for BWMs in individual pulsar timing data sets and 
for assessing the minimal projected amplitude detectable 
at a particular epoch. The timing perturbation from a 
BWM is deterministic and can be included as part of 
the timing model. The projected amplitude enters the 
timing model as a linear pa rameter and can be fit for 
using least-squares methods ( Gregory|[ 2010), yielding an 
estimate for the projected amplitude, h p , and its uncer¬ 
tainty, <Jh,p- We deal with the nonlinear parameter t B by 
searching over a grid of trial burst epochs. For all of our 
timing model fits and calculation of timing residuals, w e 
use the software package TEMP02 (Edward s et al.| 2006). 

The modification to the timing solution caused by in- 
cluding a burst of projected amplitude h p at time t B can 
be assessed by computing the likelihood ratio for a model 
with and without a burst: 


r (h p ,t B ) = exp [x 2 (h p ,t B ) -X 2 (0,£b)]^- (12) 
In this expression, 

X 2 {K^b) = R T (/ip,iB)C _1 R(/ip,t B ), (13) 

where R(h p , t B ) are the timing residuals when a burst of 
projected amplitude h p at time t B is included as part 
of the timing model. The least-squares estimator for 
the projected burst amplitude at a trial burst epoch t Bl 
h p (t B ), maximizes T at that trial burst epoch; call this 
value T (t B ). For a fixed trial burst epoch, if the data are 
consistent with the noise model, D m 21nT, the reduc¬ 
tion in the y 2 of the residuals caused by introducing a 
projected burst amplitude to the timing model fit will be 
a random variable following a y 2 distribution with one 
degree of freedom: 


hip) = (27rD)- 1 / 2 exp (—D/2). (14) 


Since the burst amplitude is a linear parameter in the 
timing model, we expect the y 2 value of the residuals to 
respond quadratically to the burst amplitude, i.e., 


X 2 (0 ^b) = X 2 ( h p (t B ),t B ) + h p (t B )/a h ^ p (t B ) 


(15) 


or D(t B ) = [h p {t B )/(Jh^ p {t B )] 2 . The cr^p values are the 

least-squares 1-cr amplitude uncertainties on e.g., the 
quantities plotted in Figure 1. 

If we compute D for two trial epochs that are very close 
together, the results will be correlated. Because of these 
correlations, when computing D along a densely sampled 
grid of N t trial burst epochs, the number of effectively 
independent trial burst epochs tested is Ni < N t . The 
probability distribution function for the maximum value 
of D along the grid, £> max , is 


fNi (Tmax) = Njfx (D max )F 1 (Anax)^”^ , (16) 
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Fig. 3. — Cumulative distributions of D max values from 1000 
simulations of noise-like timing residuals for three of the pulsars in 
our sample. Overlaid are the anticipated cumulative distributions if 
Nj , the effective number of trial burst epochs tested, is one, five, or 
twenty. For all pulsars in our sample, the analytically anticipated 
cumulative distribution best fits the results of our simulations with 
Nj between 4 and 5. 

where F\ is the cumulative distribution of /i, 

F 1 (D)=e rf(Vh/2)- (17) 

The cumulative distribution for D max associated with 
fish is then simply 

F Ni (D max) = erf"' (yi > max y . (is) 

The anticipated false alarm probability for noise alone to 
exceed a threshold value of D max , Afresh, is just 

A, (Ahresh) = 1 AAAhresh). (19) 

For a fixed allowable false alarm probability, Ahresh 
grows logarithmically with A if A 1. 

To estimate A, for each pulsar, we generated 1000 
simulated sets of TOAs that matched the real data set in 
number of TOAs, observing schedule, and timing model, 
but yielded timing residuals consistent with our noise 
models. We then fit for h p (tB ) along an equispaced grid 
of twenty trial burst epochs between MJD 53500 and 
54900 and computed l)(A) to get D max - 
In Figure 3, we show the cumulative distribution of 
the 1000 Umax values from our simulations for three pul¬ 
sars: J0030+0451 and J2145—0750, observed by Arecibo 
and the GBT, respectively, with rms timing residuals of 
~ 100 ns, and J1713+0747, observed by both Arecibo 
and the GBT with an rms timing residual of ^ 30 ns. 
We also plot the theoretically anticipated curves for A 
equal to one, five, and twenty. In fitting the anticipated 
cumulative distributions of l) max with A as a free pa¬ 
rameter to the results of our simulations, we find that for 
all pulsars, A is between four and five. As A increases, 
large values of Umax occur more frequently in the pres¬ 
ence of pure noise. We take the conservative approach 
and assume that there are five independent trial burst 


epochs to test. With A fixed at five, inverting Equa¬ 
tion 19 allows us to compute Ahresh for any desired al¬ 
lowable false alarm probability; false alarm probabilities 
of approximately five and one percent are expected for 
Ahresh equal to 6.60 and 9.54, respectively. 

6. SEARCHES FOR BWMS IN THE EARTH TERM 

Searches for BWMs in the Earth term have several ad¬ 
vantages over searches in pulsar terms. The timing per¬ 
turbation from a BWM will turn on simultaneously for all 
pulsars in the PTA if it occurs in the Earth term, provid¬ 
ing a powerful means by which pulsar-specific phenomena 
(e.g., glitches) can be ruled out. Variations in the pro¬ 
jected BWM amplitude from pulsar to pulsar provide 
information about the location of the GW source and its 
polarization. As such, for a trial BWM source location, 
polarization, and epoch, the residuals of all pulsars in the 
PTA can be combined into a coherent global fit for the 
true burst amplitude , A- Such global fits, described in 
Madison et al.| (|2014|), have better amplitude sensitivity 
than what is attainable with any one timing data set and 
can be carried out with a similar least-squares apparatus 
as is used in pulsar term searches. 

Global least-squares fitting techniques were rece ntly 
used by the PPTA to search for BWMs (|Wang et al.' 
2015). The PPTA data set those authors analyzed con- 


tamed many fewer TOAs than the NANOGrav data set 
we are considering here owing to NANOGrav’s practice 
of reporting many TOAs from a single observing epoch 
but from different frequency channels. With these multi¬ 
frequency TOAs, NANOGrav includes as part of its tim¬ 
ing models numerous chromatic parameters to account 
for profile evolution across our wide-bandwidth receivers 
and time-variable DM, so, we fit for many more timing 
model parameters than the PPTA. 

Suppose N is the number of timing model parame¬ 
ters being fit for each of M pulsars in a PTA (in reality, 
N varies from pulsar to pulsar). Then carrying out a 
global fit for the amplitude of a BWM requires the in¬ 
version of an (NM + 1) x (NM + 1) matrix, a proce¬ 
dure that requires computation time oc (iVM + 1) 3 . This 
global fit has to be done over AAA trials in the four¬ 
dimensional phase space of sky position, burst epoch, and 
burst polarization (At A,, an d A are, respectively, the 
number of sky positions, polarization angles, and burst 
epochs tested). Methods requiring these global fits are 
thus computationally onerous for NANOGrav because of 
the large number of timing model parameters for which 
we fit. More importantly, the total number of suitably 
stable MSPs being timed by NANOGrav is growing with 
time as projects like the Arecibo PALFA survey (e.g. 


Cordes et al.| 2006 Swiggum et al. 2014) , the Arecibo 
All-sk y 327-MHz Drift Pulsar Survey (e.g., Deneva et al. 
2013|), an d the GBT Norther n Celestial Cap Pulsar Sur¬ 
vey (e.g., Stovall et al.||2014 ) continue to find more pul¬ 


sars. A search procedure based on global fits will not 
scale well to future PTA endeavors. We use a different 
technique that avoids having to carry out any global fits 
and is computationally much more efficient. 

6.1. Accelerated Earth-term Searches 

Over a five-dimensional grid of trial burst sky posi¬ 
tions, Qij, polarizations, epochs, A,z, and ampli- 
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tudes, we search for the maximum in the global 

likelihood ratio surface, 


r hB,m) (20) 

M 

K= 1 

where Bk is the projection factor, £>(#,</>), for the Kth 
pulsar. The Tk quantities are analogous to the likelihood 
ratios in Equation 12. The angles 0 and </> depend on the 
location of the Kth pulsar and the trial burst location, 
Qij. The trial burst polarization angle, -0^, influences </>. 
The total number of trials in this five-dimensional search 
is NnN^N t Nh, where Nh is the number of trial burst 
amplitudes tested. 

Ostensibly, the procedure appears to require that for 
each grid point in the five-dimensional space and for 
each pulsar in our sample, we incorporate the signa¬ 
ture of a BWM of projected amplitude Bx^ij, 
occurring at time tsg into the timing model of the 
Kth pulsar, refit the timing model, and compute T^- 
This would take computation time oc NnN^N t NhMN 3 , 
a speed up over the global fitting scheme so long as 
Nh < (NM + 1) 3 /MN 3 . The average N for our sam¬ 
ple of M = 12 pulsars is approximately 90, meaning Nh 
needs to be less than approximately 140 for this tech¬ 
nique to be faster. 

However, there is a greater speed-up to be had. Varia¬ 
tions in Slij, 'ipk, and only alter the projected BWM 
amplitude along different pulsar lines of sight. We can 
precompute the two-dimensional (projected burst ampli¬ 
tude and burst epoch) likelihood ratio surface for each 
pulsar in our sample, a procedure requiring computa¬ 
tion time oc N p N t MN 3 , where N p is the number of trial 
projected burst amplitudes considered. We then con¬ 
struct the full five-dimensional likelihood ratio surface 
from Equation 20 by figuring out what the projected 
burst amplitude in the Kth pulsar would be for a partic¬ 
ular choice of fbj, t/Tc, and looking to the pre¬ 

computed two-dimensional likelihood ratio surface for 
the Kth pulsar, and interpolating between the nearest 
projected burst amplitudes tested and the projected am¬ 
plitude of interest from the five-dimensional search; this 
step requires no additional timing model fits (the com¬ 
putationally costly step), is very fast, and is essentially 
an exercise in efficiently searching lookup tables. 

Using precomputed two-dimensional likelihood ratio 
surfaces for a global Earth-term BWM search as we have 
just described will lead to a speed up over the global 
fitting scheme if N p < NqN^(NM + 1) 3 /MJV 3 . For 
our search, we have tested N t = 40 trial burst epochs 
evenly spaced between MJDs 53541 and 54995, N.^ = 17 
trial burst polarization angles evenly spaced between 0 
and 7r, and Nq = 1598 trial burst sky positions isotrop¬ 
ically distributed on the sky. As a brief but important 
aside, we have chosen to search within the particular win¬ 
dow of dates just mentioned because for each pulsar in 
our sample, there is at least one collection of multifre¬ 
quency observations before and after this window, with 
the limiting pulsar on the early side being J0613—0200 
and J2145—0750 on the late side. With the grid pa¬ 
rameters we have chosen, so long as N p is less than ap¬ 
proximately 4 x 10 6 , using precomputed two-dimensional 


likelihood ratio surfaces will be faster than global fit¬ 
ting schemes. In practice, we have set N p = 300, test¬ 
ing 150 trial projected amplitudes logarithmically spaced 
between 5 x 10“ 17 and 10“ 12 and their negatives. Uti¬ 
lizing the precomputed two-dimensional likelihood ratio 
surfaces is thus faster than using global fitting by a fac¬ 
tor of over thirteen thousand and is the only reason why 
searching such a densely sampled grid is feasible. 

Once the global, five-dimensional likelihood ratio sur¬ 
face is computed, for fixed and ts,i, we isolate 

the value of that maximizes T^; we call it hs- This 
four-dimensional surface, which we call is approxi¬ 
mately equal to Yq , what we would get if we had carried 
out the computationally costly global fits we have taken 
great care to avoid, differing from Tq only because of 
the finite resolution of our trial burst amplitude grid. 
We can then, as in the pulsar-term case, consider the 
false alarm statistics of the quantity Dq — 2 In T^. For 
any fixed grid point in our four-dimensional parameter 
space, if the data are consistent with our noise models, 
Dq again follows y 2 statistics with one degree of free¬ 
dom. With noise-like data, if we consider the probabil¬ 
ity distribution of Dq, max, the maximum value of Dq 
over the whole grid, also follows Equation 16, but with 
a larger number of effective independent trials sampled 
than in the pulsar-term case owing to a search not just 
over time, but also over polarization angle and sky posi¬ 
tion. We call the effective number of independent trials 
in the global search Nq. We discuss Nq in more detail 
in Section 8. 


6.2. Assessing BWM Amplitude Uncertainty 

A similar relationship exists between Dg, h#, and ah 
as exists between £), h p and ah :P as discussed in Section 

5, i.e., D = (h p /ah , p ) 2 . Since T^ is just the product 
of all the pulsar-specific likelihood ratio surfaces (with 
appropriate projection factors, as in Equation 20), 


M 


M 


Dg = D K {B K hB ) = h 2 B ^(B K /*h, P ,K) 2 . (21) 

K =1 K=l 

So, we have 

Dg = (h B /(Th) 2 i ( 22 ) 

where 


°h = 


M 


-.- 1/2 


y ~2{BK/vh, P , K ) 2 


,K =1 


(23) 


Equ ation 2 3 is eq u ivalen t to Equation 32 from |van| 
Haasteren & Levin (2010). While those authors ana¬ 
lyzed idealized pulsar timing data sets with equispaced 
observing epochs, uniform TO A uncertainties, and very 
simple timing models, their result holds more generally 
and applies here. 

6.3. Cross-checks with Bayesian Methods 

As an independe nt check, we also carr y out the 
Bayesian method of van Haasteren & Levin (2010), im¬ 
plemented in the software package Piccarcf 27 ] developed 

21 https://github.com/vhaasteren/piccard 










TABLE 1 

Summary of Results from Pulsar Term BWM Search 


r 


PSR 

^max/iQ - 13 

./.max 

l B,i 

(MJD) 

D max 

JoglO Us) 

J0030+0451 

-0.87T0.77 

54100 

1.2 

-0.11 

J0613—0200 

0.68T0.50 

54584 

1.8 

-0.20 

J1012+5307 

1.47T0.79 

54062 

3.4 

-0.54 

J1455—3330 

-2.57i3.35 

54211 

0.5 

-0.02 

J1640+2224 

-2.65il.45 

53801 

3.3 

-0.52 

J1713+0747 

0.08T0.07 

54808 

1.4 

-0.13 

J1744—1134 

—2.06=11.19 

53541 

2.9 

-0.45 

B1855+09 

1.00i0.73 

53578 

1.8 

-0.21 

J1909—3744 

— 1.29=10.74 

54994 

3.0 

-0.45 

J1918—0642 

-8.29T3.39 

54994 

5.9 

-1.15 

J2145—0750 

26.60=111.51 

54994 

5.3 

-0.99 

J2317+1439 

—2.404=1.20 

53541 

3.9 

-0.67 


Note. — Description of most significant BWM-like signal de¬ 
tected in the pulsar term of each of 12 NANOGrav data sets. From 
left to right, the columns are the name of the pulsar, the most sig¬ 
nificant projected BWM amplitude calculated from least-squares 
fitting, the trial burst epoch at which the most significant ampli¬ 
tude was found, the corresponding Z) max value, and the logarithm 
of the anticipated false alarm probability for that value of D max 
assuming JVj, the number of effectively independent trial burst 
epochs tested, is five. 


independently of the Madison et al. (2014) method. In 
brief, this Bayesian method takes the likelihood of Equa¬ 
tion 9 as a function of the noise parameters E, the tim¬ 
ing model parameters Jp, and the BWM parameters 
(Qij 5 'i/jk 5 t]3,i i hB,m) • Similar to what we do for our fre¬ 
quent ist search, we ke ep the noise p arame ters fixed to 


the values obtained in Arzoumanian et al. (2014). Our 
prior distributions are hat in all stated BWM parame¬ 
ters, where we note that the prior in is taken flat 
over the sphere. We have been able to confirm that the 
results we describe in this paper are identical for the fre- 
quentist and Bayesian methods. 


7. PULSAR TERM RESULTS 

For all twelve of the pulsar data sets we analyze, our 
search for pulsar-term BWMs occurring between MJDs 
53541 and 54994 yields results that are entirely consis¬ 
tent with our noise models. In Table 1, we summarize 
the key results in our search for pulsar-term BWMs. For 
each pulsar in our sample, we list the pulsar name, the 
most significant value of h p (according to the D num¬ 
ber) along with its 1-cr uncertainty, the trial burst epoch 
at which this most significant amplitude was found, the 
Dmax value for that pulsar, and the logarithm of the 
false alarm probability anticipated from noise alone for 
that value of l) max assuming Nj = 5. The most sig¬ 
nificant event we find is in the data for J1918—0642; it 
is consistent with approximately seven percent of noise 
realizations. 

The epochs at which the most significant projected 
BWM amplitudes occur are distributed nearly uniformly 
throughout the range of trial epochs we tested with re¬ 
peat values occurring only at the very first and last 
epochs tested. Clustering at the edges of the window 
of tested epochs is not surprising. If there are few tim¬ 
ing residuals outside of the window of tested trial burst 
epochs, when testing the first or last trial epoch, the 
pre- or post-burst timing model is constrained by a small 
number of data points. If the residuals at the edge of the 
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Fig. 4.— Results of a search for a BWM in the NANOGrav 
five-year data set for PSR J1713+0747, the single most sensitive 
data set in our sample for such searches. The bottom panel shows 
the epoch-averaged timing residuals for J1713+0747 as they are in 
Figure 2. The second panel from the bottom shows for the 

40 trial burst epochs we tested. The third panel from the bottom 
shows the detection statistic D varying over the span of trial burst 
epochs tested. The remaining panels show a “heat map” of the 
two-dimensional likelihood ratio surface, as in Equation 12, that 
we use for our Earth-term analysis. 


data set are slight outliers, allowing for an instantaneous 
change in the spin period of the pulsar near the begin¬ 
ning or end of the data set can bring them more in line 
with zero and lead to a modest reduction of the y 2 value 
of the residuals. 

The individual data set in our sample that most tightly 
constrains BWMs is for J1713+0747. Th is is in line with 
the expectations of Madison et al. (2014). In Figure 4, we 
show the detailed BWM search results for J1713+0747. 
The bottom panel is simply the epoch-averaged residuals 
as in Figure 2; we plot them again here to explicitly show 
how the interval of trial burst epochs tested overlaps with 
the TOA coverage. The second panel from the bottom 
shows hp ± (Jh,p at the 40 trial burst epochs we tested. 
The best-fit projected amplitude is completely consistent 
with zero at each trial burst epoch considered. Further 
echoing this, in the third panel from the bottom, we show 
the value of D as it varies with trial burst epoch. Never 
does D approach the values necessary to be inconsistent 
with noise at the five- or one-percent level. In the re¬ 
maining panels of Figure 4, we show a “heat map” of 
the two-dimensional likelihood ratio surface, as in Equa¬ 
tion 12, that we use in our Earth-term search. 






































































9 


Epoch = 54659 Pol. Angle = 1.37 





Fig. 5.— A time and polarization slice of our Earth-term search. 
The plotted quantity, Dq , is equivalent to the best-fit BWM 

amplitude, Hb, in units of the 1-cr uncertainty on the amplitude, 
cTh . At the location indicated by the black circle, cr^ = 6.07x10“ 15 , 
the smallest value of in our entire search. In our whole search, 

the maximum value of Dq we find is approximately 2.46, entirely 
consistent with our noise models at the 95% level even if Nq is 
no greater than 5 as we used in our pulsar-term searches. The 
diamonds indicate the positions of the 12 pulsars in our analysis. 
The largest diamond represents J1713+0747. 



RA 


Fig. 6.— The maximum luminosity distance of a SMBHB merger 
causing a BWM detectable with 95% confidence given our sensi¬ 
tivity averaged over burst epoch and polarization angle. We have 
assumed the binary had a typical inclination angle of 7r/3 and a 
reduced mass of 10 9 M©; the plotted distances scale linearly with 
this fiducial reduced mass. The diamonds indicate the positions 
of the 12 pulsars in our analysis. The largest diamond represents 
J1713+0747. The green triangle represents the position of the cen¬ 
ter of the Virgo Cluster. Just 16.5 Mpc from Earth, the Virgo 
Cluster is near the edge of the volume in which we can detect 
BWMs with these properties. 

8. EARTH TERM RESULTS 

In Figure 5, we show a two-dimensional slice of our 
Earth-term search results. For fixed trial burst polar- 

ization angle and epoch, we depict the quantity Dq as 
it varies with trial burst source location. The quantity 

Dq 2 is the best-fit BWM amplitude, h#, in units of the 
uncertainty on the amplitude, ah (see Equation 22). The 
slice we show contains the smallest value of ah found in 
our entire search, i.e., the point in our parameter space 
where we are most sensitive to a BWM; its location on 
the sky is indicated by the black circle very near the lo¬ 


cation of J1713+0747. 

The single largest value of Dq we find is A^max = 
6.03, which corresponds to Hq = 2.46(7^. We mentioned 
at the end of Section 6.1 that the effective number of 
trials tested in an Earth-term BWM search, Nq , will 
be larger than Nj = 5 because of the search over not 
just many trial burst epochs, but also over many trial 
burst source locations and polarizations, implying that 
for a fixed allowable false-alarm probability we have to 
raise the threshold we impose on Dq above the threshold 
used on the detection statistic in a pulsar-term search. 
However, D^max is small enough that even if Nq were 
only five, this result would be entirely consistent with 
more than 95% of realizations of our noise. _ 

In the BWM analysis carried out by Wang et al. (2015), 
when faced with marginally high values of their test 
statistic, they conducted an extensive suite of simula¬ 
tions to assess Nq, comparable to what we have done 
to assess Nj (and depicted in Figure 3), and were able 
to justifiably raise the detection threshold on their test 
statistic and rule out a detection. We find no comparably 
large values of Dq that exceed our detection threshold 
even if we underestimate Nq as five. 

We do still want an estimate of Nq as it will allow us to 
apply upper li mits to the popu l atio n of B WMs given our 
non-detection. Cornish & van Haasteren (2014) recently 
demonstrated that the response of a ETA to GWs of any 
waveform can be decomposed into a linear combination 
of a finite number of modes, or sky maps. The number 
of modes required is equal to twice the number of pulsars 
in the array (the factor of two accounts for the two pos¬ 
sible polarization modes of GWs). We use their result 
to estimate that the number of statistically independent 
samples in our Earth-term search over source-position 
and polarization space is 24, twice the number of pulsars 
used in our analysis. Given five independent samples in 
time, we thus adopt Nq = 120. 

Adopting Nq = 120 is a conservative estimate as our 
sensitivity to BWMs is so strongly dominated by a sin¬ 
gle pulsar, so the effective number of pulsars in our array 
is fewer than 12. As mentioned following Equation 19, 
for a fixed allowable false alarm probability, if Nq 1, 
Dq ,thresh only diverges logarithmically with Nq , so over¬ 
estimates of Nq are not exceedingly deleterious for the 
purpose of setting upper limits. With Nq = 120, setting 
Dq : thresh = 12.4 assures a false alarm probability of less 
than five percent. This is equivalent to requiring that Hb 
be greater than 3.52(7^ in order for it to be inconsistent 
with 95% of realizations of noise. Again, in our Earth- 
term search, we find no signal that meets or exceeds this 
level of significance. 

In Figure 6, we have averaged ah over trial burst po¬ 
larization angles and epochs and shown the maximum 
luminosity distance at which a SMBHB merger with an 
inclination angle of 7r/3 and /i = 10 9 M© (consistent with 
Equation 2) would be detectable with our data set with 
95% confidence, or where Jib = 3.52(cr/ l )^ ? t- Our sen¬ 
sitivity is worst near the position of PSR J0613—0200. 
This has little to do with J0613—0200, but is instead 
because these sky positions are antipodal to our great¬ 
est concentration of pulsars, especially our most sensitive 
pulsar, J1713+0747. The green triangle in Figure 6 in- 















10 


dicates the position of the Virgo Cluster. Just 16.5 Mpc 
from Earth, the Virgo Cluster falls very near the edge of 
the volume over which we are sensitive to BWMs from 
/i = 10 9 M 0 , X = 7r/3 binary mergers. 

9. BWM RATE-AMPLITUDE CONSTRAINTS 


To synthesize the results from both our Earth- and 
pulsar-term analyses, we derive constraints on A(> ft#), 
the annual rate of BWMs from any part of the sky with 
any polarization having amplitudes greater than he, as¬ 
suming they occur as a Poisson process. This quantity 
is readily relatable to astrophysical mod els of p rocesses 
produ cing BWMs, i.e ., Equation 15 of|Cordes fc Jenet] 
(|2012 |), Equation 21 of Madison et al. (2014), or Equation 
17 of | Wang et al. (2015|). 

Toward this end, define the quantities, 


TE(hB,n )= 


AtAQA'iJj 

47r 2 


N t N n 

EEE Q(Hb — n(T hj ijk), 

i j k 


(24) 


and 


rp{h B ,n) = ^ E E/ e ) d * da ’ 

(25) 

where At, Af}, and Aip describe the grid spacing in the 
Earth-term search. The quantities te and rp are the 
total time that the PTA had n-a sensitivity to a BWM of 
amplitude at least hp in the Earth-term and pulsar-term, 
respectively, weighted by the fraction of the total source- 
location and polarization angle space over which that 
sensitivity was achieved. Our definitions for te an d rp 
differ from nearly identical definitions in Mad ison et al.| 
(2014|) by an overall factor of two. In that work, it was 
mistakenly assumed that only BWM polarization angles 
between 0 and 7 t/ 2 must be considered. We here correctly 
carry out a search that tests polarization angles between 
0 and 7r. 

With the definitions for rp and rp in place, we can 
derive constraints on A(> hp) from our Earth-term and 
pulsar-term analyses: 


M N t 


A (>M<~ ln(1 (26) 

Ti 

where Q is the probability that at least one BWM oc¬ 
curring at rate A(> hp) encounters the PTA during the 
time Ti and i is a placeholder for either E or P. 

In Figure 7, we show our constraints on A(> hp) 
from our Earth- and pulsar-term analyses. We have set 
Q = 0.95. We have set n = 3.52 for both our Earth- and 
pulsar-term constraints. This number comes from requir¬ 
ing a false alarm probability of less than five-percent in 
an Earth-term search with Nq = 120. Our estimate of 
Nq is likely significantly larger than it needs to be. A 
detailed suite of simulations could give us a more realistic 
assessment of Nq which would lead to more constrain¬ 
ing upper limits, but for a fixed allowable false alarm 
probability, the amplitude a BWM must exceed to be 
inconsistent with the noise scales as the square root of 
the natural logarithm of Nq, so the improvement to the 
upper limit from this analysis would be very slight. 



Fig. 7.— 95% confidence Earth-term upper bound on the rate, 
A(> Hb), of BWMs occurring with amplitudes at or above ampli¬ 
tudes Hb- The two curves come from our Earth-term and pulsar- 
term analyses. The pulsar-term probes lower-rate events at high 
amplitudes because individual pulsar terms contain causally inde¬ 
pendent information. The near total convergence of the Earth- and 
pulsar-term curves at low BWM amplitudes demonstrates that a 
single pulsar, J1713+0747, is dominating our sensitivity. 


Since there are fewer effectively independent trials in a 
pulsar-term search, we should use a lower value of n for a 
95% confidence upper limit. But, by treating the Earth- 
and pulsar-terms comparably, we are able to highlight 
an important fact about our data set. The near-total 
convergence of our Earth- and pulsar-term constraints 
for low values of hp indicates that our upper limits are 
almost entirely dominated by a single pul sar data set: 
J1713+0747. This was anticipated by |Madison et al. 


(2014) and holds true even with the more sophisticated 
noise modeling we have employed in this analysis. 


10. CONCLUSION 

In this paper, we have conducted a search for BWMs 
in the first five years of NANOGrav data. We did not 
detect any BWMs. Based on our current understanding 
of SMBHBs, the most conventional anticipated so urce of 
bright BWM s, it is unsurprising that we did not. |Wang| 
et al. (2015) predict that BWMs from SMBHB mergers 
exceeding amplitudes of 10 ~ 14 occur at a rate of i nst a 


few every 10 5 yr. However, Cordes & Jenet (2012) con¬ 
clude that because of large uncertainties in things such as 
the inspiral rate of SMBHBs, the actual rate of BWMs is 
essentially unconstrained. Also, we stress that memory 
is a very general feature of bright GW events and large- 
amplitude BWMs may be produced by wholly unex¬ 
pected phenomena occurring at unconstrained rates; on¬ 
going searches for BWMs are crucial for utilizing t he raw 


d isco very potential of PTA observations (Cutler et al. 
2014). The methods developed in this paper are reacf 


ily generalizable to future, more sensitive data sets and 
the accelerated search techniques we have developed for 
Earth-term analyses will greatly expedite future BWM 
searches. 


Our Figure 7 and Figure 10 from Wang et al. (2015) 
show constraints on the rate of BVVlVls from initial 
NANOGrav and PPTA data releases, respectively. The 
NANOGrav constraints and the PPTA constraints are 
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quite similar. Both probe to BWM amplitudes of ap¬ 
proximately 2 x 10 -14 ; the PPTA upper limit extends 
to lower rates than the NANOGrav upper limit in part 
because the PPTA analyzed a slightly longer data span 
than we have a nalyzed here. Our Figure 6 and Fig¬ 
ure 9 from Wang et al. (2015) illustrate the sensi¬ 
tivity of NANOGrav and the PPTA to BWMs as it 
varies over the sky. Combined, these figures help sup¬ 
port the science case of the IPTA. The area of the 
sky where NANOGrav’s sensitivity to BWMs is worst 
is non-concentric with the area of the sky where the 
PPTA’s sensitivity is worst. Joint analysis of data from 
NANOGrav and the PPTA will lead to more uniform 
sensitivity to BWMs over the whole sky and more con¬ 
straining upper limits. Though the EPTA has not yet 
conducted a search for BWMs with their data, the inclu¬ 
sion of their data in the joint IPTA data set will likely 
play a similarly important role in improving the unifor¬ 
mity of the IPTA’s BWM sensitivity. 

NANOGrav has made great strides in improving its 
BWM sensitivity since the collection of this initial data 
set. The pulsar timing backends at Arecibo and Green 
Bank, ASP and GASP, have been upgraded to the Puerto 
Rican Ultimate Pulsar Processing Instrument (PUPPI) 
and the Green Bank Ultimate Pulsa r Processing Instru¬ 


ment (GUPPI; DuPlain et al. 
ceptionally wide bandwidths 


2008), backends with ex- 
ratThave reduced the rms 


timing errors on most pulsars being timed by a factor 
of two to three. Furthermore, NANOGrav is now reg¬ 
ularly timing more than 40 pulsars rather than just 17 
as in the first five years. Finally, just having a longer 
timing baseline on the pulsars we have analyzed in this 
paper is a great boon to our BWM sensitivity. All else 
being equal, sensitivity to BWMs scales approximately 


as T -3 / 2 , where T is the span of the data set. If we then 
focus on SMBHBs above a minimum reduced mass (as 
in Figure 6), the volume of space in which we are sensi¬ 
tive to memory from their mergers grows as T 9 / 2 ; as the 
volume in which we are sensitive to billion solar mass 
mergers already encompasses parts of the Virgo cluster, 
this strong scaling of volume probed with time means 


that many more astrophysically interesting systems will 
enter our detection horizon with continued PTA observa¬ 
tions. NANOGrav is preparing a collection of nine years 
worth of data that will be a significantly more sensitive 
probe of BWMs than any PTA data set that has yet been 
analyzed. We plan to apply the techniques used in this 
paper to the nine-year NANOGrav data set to produce 
unprecedented constraints on BWMs. 
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APPENDIX 


TABLE 2 

Table of Key Symbols 


Symbol 

Description 

Dimensions 

■ 

a “hat” on any quantity indicates that the BWM amplitude used to calculate it is 
the result of least-squares fitting 

dimensionless 


a “tilde” on any quantity indicates that the BWM amplitude used to calculate it maximizes 
the global likelihood (i.e., Equation 20) among trial amplitudes in a grid search 

dimensionless 

b(0,4>) ■ 

projection factor accounting for the geometric configuration of the Earth, 
pulsar, and burst source and the burst polarization angle 

dimensionless 

(-* , , 

,rr' 

covariance from noise between a TO A from epoch r and frequency 
channel v and a TOA from epoch t' and frequency channel v' 

time 2 

D(h p , t B ) 

21nT (h p ,t B ) 

dimensionless 

d g 

• • 2 In r G 

dimensionless 

d l 

luminosity distance from Earth to BWM source 

distance 

E 

EFAC, a constant multiplier on the rms perturbation attributed to radiometer noise 

dimensionless 

fk 

the p.d.f. for the maximum of k random numbers drawn from a % 2 
distribution with one degree of freedom 

dimensionless 

F k 

the c.d.f. associated with /& 

dimensionless 


one minus the c.d.f. associated with 

dimensionless 

hs 

amplitude of a BWM 

dimensionless 

hp 

projected BWM amplitude, 

dimensionless 

J 

jitter-like covariance between TO As from the same epoch but different frequency channels 

time 
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TABLE 2 — Continued 


Symbol 


Description 


Dimensions 


N : 

N g 

Q 

to 
11 
tB 

u 

T(h p , ts) 

r G 

At 

A(> Hb) 

°h 

&h,p 

°J 

^S/N 

T E (h B ,n) 

rp(h B ,n) 

X 2 (h P ,t B ) 


the effective number of trial BWM epochs tested in a single-pulsar BWM search 
the effective number of trials in the space of BWM epoch, sky position, 
and polarization angle in a global Earth-term BWM search 

EQUAD, the rms of a Gaussian noise process added in quadrature to radiometer noise 
time at which BWM wavefront encounters the Earth 

time at which BWM wavefront is observed from Earth to encounter a pulsar 

epoch BWM is observed to occur, either to or t\ 

the “exploder” matrix that map s observatio n epochs to the full set 

of TOAs, as discussed in Arzoumanian et al. |2014) 

likelihood of a pulsar’s TOAs assuming a BWM of projected amplitude h p occurred 
at epoch ts divided by the likelihood of the TOAs assuming no BWM occurred 
the global likelihood ratio, or product of pulsar-wise likelihood ratios (i.e., Equation 20) 
for a trial BWM of fixed sky location, polarization, epoch, and amplitude 
perturbation to pulse times of arrival 

the rate that BWMs with amplitudes greater than hs encounter our PTA 

reduced mass of a SMBHB, Mi M 2 /(Mi + M 2 ) 

uncertainty on the amplitude of a BWM 

uncertainty on the projected amplitude of a BWM 

rms timing perturbation from pulse phase jitter noise 

rms timing perturbation from radiometer noise 

total time that an Earth-term search could yield n-cr^ detections of 
a BWM with an amplitude greater than hs 

total time that a pulsar-term search could yield n-cr^ detections of 
a BWM with an amplitude greater than hs 

the square of the timing residuals when a BWM is included in the timing model 
weighted by the inverse noise covariance matrix (see Equation 13) 


dimensionless 

dimensionless 

time 

time 

time 

time 

dimensionless 

dimensionless 

dimensionless 

time 

time -1 

mass 

dimensionless 

dimensionless 

time 

time 

time 

time 

dimensionless 


Note. — For reference, descriptions of important or frequently used symbols. 
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